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Eukaryotic cells possess motility mechanisms allowing them not only to self-propel but also to exert forces on 
obstacles (to push) and to carry cargoes (to pull). To study the inherent asymmetry between active pushing and 
pulling we model a crawling acto-myosin cell extract as a ID layer of active gel subjected to external forces. We 
show that pushing is controlled by protrusion and that the macroscopic signature of the protrusion dominated 
motility mechanism is concavity of the force velocity relation. Instead, pulling is driven by protrusion only at 
small values of the pulling force and it is replaced by contraction when the pulling force is sufficiently large. 
This leads to more complex convex-concave structure of the force velocity relation, in particular, competition 
between protrusion and contraction can produce negative mobility in a biologically relevant range. The model 
illustrates active readjustment of the force generating machinery in response to changes in the dipole structure 
of external forces. The possibility of switching between complementary active mechanisms implies that if 
necessary 'pushers' can replace 'pullers' and visa versa. 



I. INTRODUCTION 

Most of eukaryotic cells, including, for instance, fish kera- 
tocytes, self-propel by advancing the front and retracting the 
rear. A prototypical scheme of such motility includes poly- 
merization of actin, facilitated by dynamic assembly of fo- 
cal adhesions, motor-driven contraction of acto-myosin cy- 
toskeletal network, and, finally, the detachment of adhesive 
contacts followed by depolymerization of actin that closes 
the treadmilling cycle |1-7|. All three main components of 
the motility mechanism (polymerization, contraction and ad- 
hesion) are active and require intricate regulation as well as 
a continuous supply of energy. While the general crawling 
scheme described above is compatible with both oscillatory 
and steady translocation of the cell body, in this paper we fo- 
cus on steady motility modes. 

The molecular and biochemical basis of cell motility is ba- 
sically known, however, the qualitative understanding of the 
mechanical interplay between different active components is 
hidden behind complex computational schemes involved in 
modeling of cell motility |8- 14 |. In particular, the relative me- 
chanical role of contraction and protrusion in exerting forces 
on obstacles (pushing) and carrying cargoes (pulling) is usu- 
ally obscured by geometrical and chemical complexity of the 
comprehensive mathematical models. 

Protrusion is known to be the main mechanism of pushing 
which, for instance, plays dominant role in Listeria propul- 
sion |15|. Instead, contraction is believed to be crucial for 
the ability of cells to pull organelles. An inherent functional 
disparity between protrusion-contraction components of the 
motility mechanism suggests a fundamental difference in the 
structure of the force-velocity relations associated with push- 
ing and pulling. In experimental studies pushing and pulling 
are often difficult to distinguish and most of the measured 
force-velocity data are attributed to pushing fTF-TOl. 

To separate contributions of protrusion and contraction we 
use the simplest model of an active gel and view lamel- 
lipodium as a ID fluid body f20l. The actomyosin cell ex- 
tract represented by such gel is assumed to be limited by free 
boundaries where the external loads are applied (see Fig|T]l. 
Actin treadmilling also takes place on these boundaries and is 



modelled as an influx of mass at the front and its disappear- 
ance at the rear. Active contraction is represented by a spa- 
tially homogeneous prestress generated at the microscale by 
molecular motors. Adhesion is assumed to be passive and is 
modeled in this minimal setting by viscous friction on a rigid 
background. 

Our main result is that the roles of protrusion and contrac- 
tion may by interchangeable depending on the character of the 
mechanical task performed by the cell (pushing or pulling). 
We identify an experimentally observable macroscopic sig- 
nature of the dominance of each of the two mechanisms by 
demonstrating that pushing-dominated force-velocity relation 
is concave while pulling-dominated force velocity relation 
may be convex-concave with an interval of negative mobility. 





F 








co[>tincfioii 






9- 


met iibi'niie+ cort ex 





f O — Z! 

carao nmtoi pnjf 

Figure 1 . Schematic representation of an advancing lamellopodium 
subjected to a pushing force q-f- and a pulling force g_ . 

An explicit solution of the mechanical problem shows that 
in the presence of a cargo the minimal fluid model is singular 
and we regularize it by accounting for an overall stiffness. The 
latter may be ascribed either to a membrane or to the elastic 
components of the cytoskeleton ll2TH24l . Such augmentation 
removes singularities but preserves the main qualitative pre- 
dictions of the minimal model. 

To further challenge the minimal model we study the effects 
of spatially inhomogeneous (graded) adhesion, distributed 
(bulk) depolymerization and consider the dependence of con- 
tractile stress on actin density. We show that the main qual- 
itative results regarding the convexity structure of the force 
velocity relation remain unchanged. 
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While even our augmented minimal model still underrepre- 
sents some physical effects (e.g. active adhesion, transport of 
motors, complex membrane dynamics, 3D geometry, etc. ||8]- 
[141) it allows one to go beyond force velocity relations and 
study the efficiency of cargo-pulling machinery. In particular, 
we show that a competition between protrusion and contrac- 
tion can result in a bi-modal structure of the load-efficiency re- 
lation. By using the minimal model we could also compare the 
conventional kinematic mode of driving through given poly- 
merization/depolymerization velocities with a direct control 
of external energy supplies responsible for protrusion and con- 
traction. 

Vis-a-vis the general behavior of active media, we have 
shown that an interplay between 'pushers' and 'pullers' ll25l - 
|27| can lead to observable effects in the presence of applied 
loads. The importance of the idea that different active mech- 
anisms can swap roles depending on the task goes far beyond 
the subject of cell motility. 

The paper is organized as follows. In Section II we for- 
mulate the minimal model, find explicit traveling wave solu- 
tions describing steadily advancing cells and study stability 
of these solutions. By solving the associated transport prob- 
lem we reconstruct actin density profiles in different load- 
ing regimes and reveal the mechanism behind the possibility 
of infinite density localization. We then study the distribu- 
tion of active force dipoles in the moving cell and present 
an interpretation of the negative motility regime in terms of 
a crossover between protrusion-dominated and contraction- 
dominated regimes. Elastic regularization of the minimal 
model is introduced in Section III where we consider sepa- 
rately the mean field (spring) model, the Kelvin- Voigt model 
and the Maxwell model. Other extensions of the minimal 
model incorporating inhomogeneous adhesion, distributed de- 
polymerization and density dependent contraction are ana- 
lyzed in Section IV. In Section V we explore the energetics of 
the protrusion-contraction mechanism and study the load de- 
pendence of its efficiency. The possibility of a non-kinematic 
driving of the moving cell is discussed in Section VI. The final 
Section VII contains our conclusions. 



II. THE MINIMAL MODEL 

Our starting point is the balance of forces in a ID layer of 
active gel placed on a rigid surface |20 , 28 1. While active dy- 
namics of adhesion complexes is notoriously complex |29|, a 
usual assumption made in the context of cell motility is that 
the time averaged tension generated by constantly engaging 
and disengaging focal adhesions is proportional to the veloc- 
ity of the retrograde flow lHJl [HI |M1 [30:-33|. If we neglect 
the bi-phasic effect IS [T] [34lj38ll and assume that the friction 
coefficient > is constant we obtain 

dxo = £.v, (1) 

where v{x, t) is the velocity and (j{x, t) is the stress. Here and 
throughout the paper we denote by da a partial derivative with 
respect to a. 



By using the constitutive model of an infinitely compress- 
ible viscous active fluid we can write Il20ll28l 

cr = X + Vdxv, (2) 

where > is a bulk viscosity and x > is a constant ac- 
tive pre-stress. This minimal constitutive description is clearly 
singular because the cell can be infinitely stretched or com- 
pressed and it is quite remarkable that this setting is already 
sufficient to capture the essence of active competition between 
treadmilling and contraction. 

Combining ([T]i and Q we obtain a second order differen- 
tial equation which we need to solve on a domain with/ree 
boundaries (i) and L [t) representing the front and the rear 
limits of a cell. To solve this problem on a domain W\t\\ fixed 
boundaries we need to impose two mechanical boundary con- 
ditions 

a{l±{t),t) = q±. 

These conditions introduce asymmetric loading which is the 
central concept of this paper In our notations g+ < corre- 
sponds to pushing (at the front) and > to pulling (at the 
rear). To find the unknown functions 1+ (t) and L (t) we need 
to impose two additional boundary conditions. The conven- 
tional choice is a pair of kinematical constraints Il20ll28]| 

v{i±{t),t) - i± = v±, 

where v+ > and w_ > are the polymerization and the 
depolymerization velocities, respectively. The prescribed sign 
of these velocities introduces implicit polarization of the cell 
which is necessary for initiation of motility in the absence of 
applied forces. 

If we now normalize length by y/rjJS,, time by 77/x and 
stress by x, we obtain a free boundary problem which de- 
pends on four dimensionless parameters. Two of them, v±, 
characterize internal driving and the other two, q±, describe 
external loading. It is natural, however, to work with a slightly 
different set of parameters. Thus, parameter 



prescribes polarity of the cell and, as we show later in the 
paper, gives the scale of the maximal velocity. The remaining 
kinematic parameter 

introduces the asymmetry between polymerization and de- 
polymerization and, as we show later in the paper, quantifies 
the degree of engagement of the contractile mechanism. It 
will also be convenient to define the resultant force 

Q = q--q+>,0 

which we assume to be positive and acting against the polar- 
ization direction induced by protrusion. We also introduce the 
force asymmetry factor 

q- +q+ 
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which characterizes the first moment of the external force dis- 
tribution. We notice that — 1 < e < 1 with e > correspond- 
ing to pulHng and e < - to pushing. 

The resuhing dynamic problem has a peculiar structure due 
to an implicit assumption about separation of time scales. 
More specifically, the neglect of inertia means that mechanical 
equilibrium is reached instantaneously at the time scale of the 
motion of the free boundaries (Stokes flow). The rate limiting 
factor is then kinetics of the free boundaries characterized by 
parameters v± that can be naively interpreted as describing the 
treadmilling process only. However, as we show later in the 
paper, only their sum Vm can be linked to treadmilling proper 
while their difference AV^ is a characteristic of contraction. 



A. Traveling wave solutions 

The transparency of the minimal model is due to the fact 
that our linear force balance equation with mechanical bound- 
ary conditions can be integrated in elementary functions as 
was first observed in Il20ll28l for a cell without cargo. When 
cargo is present the velocity profile can also be found explic- 
itly 

A_ cosh(L(t) - x) + ^+ cosh(Z+(t) - x) 



where 



smh{l+{t) - l^{t)) 



A± =±(l-Q(e±l)/2). 



(4) 



Knowledge of the spatial dependence and the use of kinematic 
boundary conditions allows one to obtain explicit equations 
for the functions l+{t) and L (t). Moreover, by using the total 
length L{t) = 1+ (<) — L it) we can obtain a closed dynamical 
problem 



L = Ay + (eQ - 2) tanh 



(5) 



After this equation is solved the position of the geometrical 
center of the cell G{t) = {l+{t) + l-{t))/2 can be found by 
integrating a decoupled equation with the known right hand 
side 



G = Vr,r - 



Q 



2tanh(L/2)' 



(6) 



To specify solutions of (jSj) and (j6]l we need to supply the ini- 
tial conditions L(0) and G(0) that also fix the initial velocity 
profile through ([3]l. 

In this paper we are interested in traveling wave (TW) solu- 
tions of (|5]l describing steadily translocating cells. These so- 
lutions correspond to stable critical points of ^ with L — Q 
that exist if and only if 



< Ay < 2 - eC 



(7) 



When these conditions are satisfied the length of the cell sta- 
bilizes as i — > oo at the value 



Loo = 2 tanh 



AF 



> 0. 



Alongside, the function G converges to a constant V given by 
the following force velocity relation 



Ay 



2Ay 



(8) 



Notice that the cell moves to the right against the load if 
V > Q and is dragged backwards by the load \fV<Q. The 
maximum velocity V* — Vm is achieved when there is no 
load Q = and the corresponding reference length will be 
denoted by L*^ = ^^^(g = 0). 

Since the TW regimes are stable only if 2— eQ > 0, pushing 
(e < 0) contributes to stability while pulling (e > 0) plays 
a destabilizing role. We also observe that at Ay = the 
loaded cell shrinks to a point while at Ay = 2 — eQ its length 
diverges. For singular solutions with L^o = oo which are only 
elevant in the case of pulling, the force velocity relation can 
be continuously extended by using (j6]l 



y - y™ - Q/2. 



(9) 



In Section III we show that these singular solutions of the min- 
imal model are physically meaningful and can be viewed as 
limits of the nonsingular solutions in the model with finite in- 
ternal stiffness. 

At large times we can characterize convergence of the ini- 
tial configuration to the TW profile (transient regime) by the 
formula 

\L{t)-L^\^e-"\ 

where the characteristic time of relaxation to the steady state, 

^ 2(2 - eQ) 
(2-eQ)2- Ay2' 

can be measured experimentally. After this time, which de- 
pends on both, the mechanical loading and the kinematic driv- 
ing the cell can be expected to acquire the velocity predicted 
by the steady force-velocity relation (|8]l-((9|. 



B. Force velocity relation 

The structure of the obtained force-velocity relation in the 
(y, Q) plane is illustrated in Fig|2]a,b. One can see that it is 
markedly different for e > (pulling) and e < (pushing). 
The main feature distinguishing pushing from pulling is the 
curvature of the force velocity relation which in the regular 
regimes (j8]l is given by 



Ay' 

and in the singular (pulling) regimes by 

5Q2 



0. 



One can see that the curvature is always negative in pushing 
regimes with e < which means that the corresponding force 
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Figure 2. The typical force-velocity relations in pure pulling (a) 
and pushing (b) regimes. Stress, velocity and density profiles corre- 
sponding to points A,B,C,D, a and /3 are shown in Fig|5]and Fig|4] 
Driving parameters are v- = 1.7 and w+ = 2. 



velocity relation is concave. Under pulling loads with e > 
the force velocity curve is convex for regular regimes and is 
linear for singular regimes. 

In the pushing regimes the force velocity curve is char- 
acterized by the stall force Q* = (1 - ^1 - 2eAVVm)/e 
and the maximum velocity V* = Vm, see Fig. |2] (b). The 
concavity of the force velocity relation in this case agrees 
with experiments lfT6l4T9ll . In the case of pulling, the force- 
velocity relation is convex for Q < Qc = {2 — AV) /e, where 
Loo < oo and is linear for Q > Q^, where — oo, 
see Figj2] (a). In the convex range the function V{Q) is 
non-monotone when AV < 1 and one can distinguish two 
regimes: the branch Q < Q„ = 1/e where the mobility is 
positive, V{Q) ^ V,n — Q/AV, and, as we show later in the 
paper, protrusion dominates, and the branch Qc > Q > Qn 
where the mobility is negative, V{Q) ~ eQ^/ {2AV) and the 
dominant active mechanism is contraction. Along the nega- 
tive mobility branch the cell elongates to support larger loads 
till the length diverges at a critical value Q — Qc- Beyond 
this value, we obtain configurations with infinitely separated 
boundary layers and mobility becomes again positive. 



C. Density distribution 

To interpret complex behavior of the force velocity relation 
in pulling regimes we first need to reconstruct the (actin) den- 
sity distribution inside the moving cell. The assumption of 
infinite compressibility allows one to decouple the problem of 



finding density distribution from the problem of determining 
stress and velocity profiles. 

After a 'statically determinate' mechanical problem ([TJ^Q 
is solved, the density p{x, t) can be obtained from the mass 
transport equation: 



dtp + d^{pv) = 0, 



(10) 



where the function v{x^ t) is given by (jsjl. Equation ( 10 1 must 
be supplemented by a single boundary condition 



p{l+{t),t)v+= p{l^{t),t)v. 



(11) 



which ensures that the exterior treadmilling mechanism con- 
serves the incoming mass flux 

m(t) = -p{l^(t),t)v- < 0. 

Then the total mass M — fl^j^^-l p{x,t)dx, is constant and 
all actin depolymerized at the rear is instantaneously re- 
polymerized at the front 



dM 



pil+{t),t)v+^ p{l_{t),t)v_ 



0. 




l-(P) Pi{x) i+(0) 



Figure 3. Schematic structure of a particle trajectory inside a cell 
as it approaches the steady state TW regime. Dotted lines indicate 
instantaneous treadmilling of particles from the rear boundary of the 
cell to its front boundary. 



Given an initial condition p{x,Q) = Pi{x) the mass trans- 
port problem inside the cell can be solved by method of char- 
acteristics. The initial density distribution prescribes the total 
mass M which can be absorbed into the scaling of p if we 
define dimensionless density p/pQ with po = M/ \/ri/S^. The 
distribution pi is transported in finite time along the character- 
istics from [/_ (0), ^+(0)] to the rear boundary of the cell. The 
arriving mass, characterized by the distribution p{l_{t), t), is 
then (instantaneously) transported by the treadmilling mecha- 
nism ([TTJ from the back of the cell Z_ {t) to the front of the cell 
l+{t). From there the mass is again transported by character- 
istics towards the rear boundary of the cell. This construction 
is then repeated indefinitely as we show in Fig[3] 
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To be more specific, consider, for instance, a characteristic 
curve X = 4>to{t) originating at /+(to)- The function <t>to{i) is 
a solution of the initial value problem 



This allows us to write the final expression for the steady state 
density profile in the form 



(12) 



The characteristic curve x — 4>to{t) reaches the back of 
the cell at t = ii which can be found from the condition 
0*0 (^i) = L(ti). The density evolution along the charac- 
teristic curve can now be recovered from the transportation 
condition 




Figure 4. Stress, velocity and density profiles for point a in Fig 
where e = ±1, Q = Q ,V = V* = 1.85 and for point /3 in Fig 
where e = —1, Q — Q* = 0.45, V = 0. Parameters v- — 1.7 and 
v+ = 2. 



In the traveling wave regime both density and velocity de- 
pend only on the co-moving coordinate y = x — Vt, with 
< J/ < Lao- In particular, p{x,t) = p{y). We must also 
have 

l_(t) = Vt 
l+{t)=Lao+Vt 

Now the mass balance equation can be integrated explicitly 
and we obtain 



p{y) = Hy) - V) 



du 



v{u) - V 



An internal configuration of a cell at zero load {Q — Q , 
V — V*), which is typical for both weak pushing and pulling, 
is shown in Fig|4](a). Similar profiles for stress and velocity 
have been already presented in 1201 l28l and here we comple- 
ment the picture by presenting the associated density profile. 
The density accumulation at the back of the cell is in agree- 
ment with the relative velocity distribution v{y) — V in the 
co-moving frame. One can see that this flow is globally ret- 
rograde with higher absolute value of velocity at the leading 
edge than at the traihng edge. 




Figure 5. Stress, velocity and density distribution inside a cell 
moving with the same velocity V = 0.95 in four different loading 
regimes indicated as A, B, C and D in Fig[2] Dashed line shows 
elasticity-regularized profiles corresponding to point C' in Figjs] 



piy) 



v{y) - V 



(13) 



where m is a constant mass flux and the function v{y) is given 
explicitly by (|3]l. Since the dimensionless total mass of the cell 
is equal to unity, we obtain 



(14) 




Adding loads generates a non zero mean flow in the 
co-moving frame and makes the profile steeper in pushing 
regimes and more shallow in pulling regimes. For instance, 
in Fig|4] (/3) we show the configuration coiTesponding to (se- 
vere) pushing at the stall force conditions Q — Q*. In this 
and similar regimes shown in Fig|5](A) and Fig|5](D) the rel- 
ative flow with respect to the average velocity is prograde at 
the rear and retrograde at the front. Instead, the distribution 
v(y) — V is globally retrograde slowing down at the rear of 
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the cell where the density profile has a maximum, see also 

Eol. 

A fundamentally different set of regimes, signifying, as we 
show in the next section, a transition from protrusion domi- 
nated to contraction dominated motility, is shown for the case 
of pulling in Fig|5](B) and Fig|5](C). Here the relative flow 
with respect to the mean velocity is retrograde at both rear 
and front and is prograde in the central part of the cell. The 
distribution v{y) — V continues to be globally retrograde with 
strongest flow at the back and at the front of the cell which 
slows down in the middle part of the cell. This velocity redis- 
tribution pushes the density maximum from the back towards 
the center of the cell. In Fig|5](C) we show that as the length 
of the cell diverges both stress and velocity profiles flatten ev- 
erywhere outside infinitely narrow the boundary layers. 



D. Pushers and pullers 

In this Section we discuss physical phenomena behind the 
observed differences in the structure of force velocity relations 
in the regimes of pushing and pulling. 

We begin with an observation ||20l |39l that the global force 
balance, stating that the applied force is resisted by the friction 
force. 




does not distinguish between pushing and pulling. To see the 
role of different active agents we need to consider the balance 
of couples and by referring to asymmetric cargo we imply dif- 
ferent signs of the dipole component of the distributed load. 

By multiplying the force balance equation ([TJ in the TW 
regime by ?/ — j 2 and integrating over the body of the cell 
we obtain 

<3^ - 7*- / {y ~ ^)v{y)dy ^ -j^ [ cr{y)dy. 

^ J^oo Jq ^ -t^oo Jo 

The first term in the left hand side 

is the moment of external forces. Since we assumed that Q > 
0, pulling is associated with a positive applied dipole while 
pushing - with a negative applied dipole. The second term on 
the left hand side 

= --,— [ {y~^)v{y)dy 

J-ioo Jo ^ 

represents frictional dipole which may have different signs. 
The integral on the right hand side defines the active dipole 
which can be also rewritten as 

Ta^-j^ f °° {l + dyv)dy. 

J^oo Jo 



This term can be further decomposed into the sum Ta =Tc + 
Tp where contraction component can be written in the TW 
regime as 

T, = 1 > 
and protrusion component as 

AV 

Tp = -— <0. 

The opposite signs of these two terms suggest that the under- 
lying active mechanisms are inherently different. By using the 
terminology of the theory of active suspensions ll25l |26 j we 
can interpret the protrusion term as representing distributed 
"pushers" and the contraction term as representing distributed 
"pullers". 

We observe that due to the presence of a contraction (posi- 
tive) force dipoles the rear boundary of the cell is pulled for- 
ward while the front boundary is pulled backward. As a com- 
pensation, contraction produces internal retrograde flow at the 
rear and pro-grade flow at the front. In contrast, protrusion 
(negative) force dipole pushes the rear of the cell backward 
while the front of the cell is pushed forward. This is compen- 
sated internally by retrograde flow at the front and pro-grade 
at the rear. These flows must be superimposed with the mean 
flow V ~ — Q/ioo which is associated solely with the total 
applied force and is therefore always retrograde. 

We can now identify separate contributions of pushers and 
pullers in building the internal configurations of the cell dis- 
played in Fig. |4] (a, f3) and in Fig. |5] (A-D). For instance, 
knowing that in the protrusion (pushers) dominated regime the 
velocity gradient must be negative (retrograde at the front and 
pro-grade at the rear, plus a constant) we can conclude that 
pushers dominate in the configurations shown in Fig. |4](a, /3) 
and in Fig|5](A, D). Similarly, if we consider regular pulling 
regimes exhibiting negative mobility at Qc > Q > Qn, see 
Fig|5](B), we can conclude that here pushers are challenged 
by pullers that enforce positive velocity gradient (retrograde 
at the rear and pro-grade at the front, plus a constant). The 
situation remains qualitatively similar in the singular pulling 
regimes illustrated in Fig|5](C). 

Therefore, we can identify the point in Fig|2](a) with 
a crossover from pushers dominated to pullers dominated 
regimes. This interpretation is supported by comparing the 
magnitudes of the two competing active couples. For instance, 
in the realistic case AV^ ~ 0.3 |20|, illustrated in Fig|2](a), a 
condition that the magnitude of the contraction couple is twice 
as big as the magnitude of the protrusion couple, \Tc\ ^ 2|rp|, 
gives the value of the force 2 — AV/ tanh{AV) = 0.97 which 
is quite close to the threshold Qn = 1/e = 1- Here it is im- 
portant to mention that at Q = 0, we have \Tc\ ~ \Tp\ which 
allows the cefl to eliminate frictional couple and achieve max- 
imum velocity. 

The observed crossover correlates with the transition from 
positive to negative mobility which also takes place at Qn. 
Negative mobility has been discussed previously in the con- 
text of individual |40-43| and interacting [44, 45 1 Brownian 
motors. The regimes where velocity of the crawling cell in- 
creases with an opposing pulling force at the rear have been 
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envisioned in (46] where negative mobility was attributed to 
the coupling between the velocity of retraction and the applied 
force v_ (Q) ll47l . In our model such coupling is absent which 
shows that negative mobility may also have a different origin. 

To make quantitative predictions we use the data from |20|: 
X = 10^ Fa, C = 5 X lO^^Pa ■ m.-^ ■ s, rj = 5 x lO^Pa ■ s, 
^ 2 and w_ = 1.7. This gives for the dimensional velocity 
of the unloaded cell {x/V^W* = 0.37x 10^''r7i-s~^ and for 
its dimensional length ^/r^^L*^ — 0.3 x IO^^to. This length 
scale is of the right order of magnitude while the velocity 
scale is at least an order of magnitude smaller than the values 
recorded for keratocytes and fibroblasts iri8ll48l . In the case 
of pure pushing e = — 1, we can use the area S = 10~^^m^ to 
obtain the dimensional value of the stall force xSQ* ~ InN 
which is realistic |16-r9|. Based on these estimates we con- 
clude that negative motility may be expected in the interval 
of pulling force values 1 — 1.7nN and this prediction can be 
tested experimentally. 

E. Formation of singularities 

Formula ([T4]i shows a possibility of the two types of degen- 
eracies associated with reaching the condition rh = 0. In such 
singular regimes the treadmilling flow becomes fully blocked. 

The first type of singular behavior takes place when the 
length of the cell diverges. Here we refer to the infinite spread- 
ing of a cell in pulling regimes with Q > Qc illustrated in 
Fig|5](C). As we show in the next Section this problem can be 
fixed if elastic stiffness is taken into consideration. 




1 2 3 04 a 02 04 

AV y/L. 
(a) (b) 



Figure 6. (a) Various pulling regimes in the parameter space 
{Vm, AV). In the domain > the cell moves against the load 
while if V < the cell is dragged by the load. Along the line XY 
the cell is static resisting the load (stall force conditions). The singu- 
lar regimes correspond to the lines ZX (infinite localization) and XR 
(infinite spreading), (b) Density localization along the path indicated 
in (a) by the solid line which ends with the formation of a singularity 
at point W. The loading is e = 1 and Q = 1.6. For AV^ = 0.3 the 
singularity is located at j/i / Loo — 0.4. 



The second type of degeneracy is associated with non in- 
tegrability of {v{y) — V)^^ even for cells with finite lengths. 
Such singularity can take place when there exists a point yo 
where v{yQ) — V. To illustrate the possibility of this type 
of singular behavior consider the whole set of pure pulling 



regimes shown in Fig|6](a). Notice that the line m = in the 
{V„i, AV) plane consists of two segments: XR and ZX. The 
segment XR is associated with infinite spreading of the cell as 
discussed above. Instead, along the segment ZX the length of 
the cell remains finite while the density localizes infinitely in 
a single point inside the cell. 




((.) (6) 



Figure 7. Sketch of particles trajectories as the steady state (TW) 
mass flux approaches zero, see Fig|6]for the related density profiles, 
(a) At small values of rh each particle spends considerable in a small 
region near a line x = yo + Vt. (b) When rh = all particle 
trajectories converge to the line x — yo + Vt which leads to a blow 
up singularity. 

To locate this point consider a regular (rh ^ 0) density pro- 
file with a local maximum at y ~ yi where dyp{yi) = 0. 
Such point can be found from the equation 

A_ sinh(?;i) A+ sinh(Loo - yi)- 

One can see that yi does not depend on Vm and if we lower 
the value of K,„ till the regime with rh = is reached (see the 
solid trajectory ending at W in Fig|6](a)) we obtain yo — yi- 
The associated phenomenon of infinite mass concentration at 
the point y ^ yi = yo is illustrated in Fig|6](b). Below the 
line (ZWX) the reconstruction of mass density is not possible 
because of the interpenetration of matter 

The phenomenon of infinite density localization can be also 
illustrated through the behavior of the characteristics (particle 
trajectories). In Fig|6]we show how characteristics in the TW 
regime concentrate as one approaches point W in Fig|6] (a). 
Even before reaching the regime W the particles spend con- 
siderable time around the line x ^ yo + Vt, see Fig|6] (a). 
At the point W where rh ^ the mass flow gets completely 
blocked as we show in Fig|6](b). Notice also that due to de- 
coupling of mechanical and mass transport problems in the 
minimal problem, the velocity field in such singular regimes 
remains regular Similar to the case of infinite spreading, the 
problem of infinite localization can be resolved if we take into 
consideration internal stiffness of the cell body. 

III. ELASTIC REGULARIZATION 

A natural way to regularize the minimal model is to intro- 
duce an intermediate-time stiffness of the cell. Such stiffness 
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prevents the unloaded cells from contraction-induced collapse 
and sets the rest length and it also keeps this length from di- 
verging in the case of pulling. 

Elasticity may be associated either with the cytoskeleton 
or with the cell membrane. Membrane and cortex elastic- 
ity can be modeled in a prototypical setting as a mean field 
elastic feedback provided by elastic springs linking different 
parts of the cell |23 , 24|. Visco-elastic properties of cytoskele- 
ton strongly depend on the characteristic time of the prob- 
lem ||49| ^5]1 and the corresponding corrections to the active 
gel model in the bulk of the cell are usually incorporated ei- 
ther in the framework of a short time (Maxwell) elastic model 
ESl [30l I52ti54l or a long time (Kelvin- Voigt) elastic model 



A. Mean field elasticity 

The simplest elastic regularization of the minimal model is 
through mean field coupling between the leading and trailing 
edges of a cell ll2T1424l . If this coupling is linear elastic, the 
applied loads become 



q± q± 



where k > is a dimensionless stiffness and Lq is a pre- 
scribed dimensionless reference length. The meaning of pa- 
rameter Lq is clear from the fact that for fc > 1 and Vm — 
AV = there exists a nontrivial static solution with L^o = 
Lo{l — 1/k) (preferred shape). 

In dynamics the steady state (TW) solution is now stable 
for all AV^ > and to find L^o [Q) one needs to solve 



2fc ^°°~^° )tanh(% 

Lq \ 2 



Then, the force velocity relation can be found from the rela- 
tion 



V{Q) = K,. - 



Q 



2tanh(^) 



and its k dependence is illustrated in FigjS] We observe that 
independently of the value of k all force-velocity curves cross 
at Q = where V ^ V* . The second common intersection 
point at 



2 

V* 









1 e = 1 1 




















1 







4 

Q 




Figure 8. Force velocity relations in pure pushing and pulling modes 
with different fci,2,3,4 = {0, 0.01, 0.1, 1} and Lq = 1. Driving pa- 
rameters are v- — 1.7 and v+ — 2. Internal profiles corresponding 
to points C and C' are presented in Fig|5](C). The minimal model is 
recovered at fc = 



Positive motility 



Negative motility 



Figure 9. Domain of negative mobility in the parameter space 
(fc, AV). The boundary between regimes with positive and negative 
mobility is given by the function fc = k*{AV). 



exists when e > and Lq > L*^. 

The salient feature of the regularized model is that at fc — >^ 
the mean field force velocity curves approach their minimal 
model counterparts including both the regular regimes with 
finite cell lengths and the singular regimes with infinite cell 
lengths. However, despite similarity in shape between the 
force velocity curves in the minimal model and in the regular- 
ized model with fc ^ 0, the length of the cell in the regularized 



model is always finite so that infinite stretching, undermining 
the minimal model, does not take place. Unfortunately, by 
using the mean field elasticity model one cannot also remove 
the localization singularity in the minimal model because in 
this regularized setting the force balance remains independent 
from the mass balance. 

The phenomenon of negative mobility for the pulled cells 
survives in the mean field model and disappears only at a crit- 
ical value of the stiffness fc — k*{AV), see Fig. |9] The qual- 
itative difference in convexity between pulling and pushing 
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persists beyond k*{AV), see for instance regime with fc = 1 
in Figjs] However, at A: ^ k*{AV) the force-velocity re- 
lations associated with pushing and pulling regimes become 
similar. 

The available data on static configurations ll2Tlj24ll suggest 
that the dimensionless parameter k (normalized by x) must be 
in the range 1 — 10 which apparently excludes the negative 
mobility regimes. However, many models of cell dynamics 
are built under the assumption that long time elasticity is neg- 
ligible and essentially assume that fc = Il20ll30l 1521 1541 The 
ambiguity maybe due to the ability of cytoskeleton to fluidize 
by engaging active cross-linkers that can modify its stiffness 
over at least two orders of magnitude |57-59|. 

In view of such broad rheological flexibility of the cy- 
toskeleton, the effective stiffness may easily reach below the 
threshold fc = 1 which means that negative motility regimes 
cannot be excluded in vivo and can be probably artificially 
engineered in vitro through partial suppression of the stiff- 
ening components of acto-myosin network. If we use nu- 



merical values of parameters from Section II A we find that 
xSQi — l.3nN which gives the scale of pulling forces where 
negative mobility can be expected. Notice that this value is 
above the resolution of an atomic force microscopy cantilever 
which has been previously used in the measurements of force 
velocity curves I.17J . 



B. Kelvin- Voigt elasticity 

Instead of using spring based elastic regularization consid- 
ered above we can directly incorporate distributed elasticity 
into the constitutive model. The closest to the mean field 
model is the Kelvin- Voigt model accounting for the elastic 
response at long time scales. In the ID setting we need to 
assume that 

(7 = x + vdxV-p{p), 

where p{p) is the stress-density relation. In our version of 
Kelvin- Voigt model we further assume that this relation is lin- 
ear 



Pip) = E 



where pr is the reference density and E the elastic modulus 
ll49l ISTl . The resulting system of coupled non-dimensional 
equations can be written as 



dtP + d^ipd^a) = 
-dxxO- + a- = 1 - K{p/pr ~ 1), 



where we introduced two new nondimensional parameters 
K — E/x and pr = Pr/ Po- To find the steady state (TW) 
regimes we need to solve a simpler system 



-a = l-K{l~l) 




Figure 10. Force-velocity curves for the Kelvin- Voigt model with 
different /•f 1.2, 3, 4 = {0,0.1,1,5}. Other parameters are pr = 1, 
u_ = 1.7 and v+ = 2. The minimal model is recovered ai K = Q. 



with mechanical boundary conditions, f7(0) = g_ and 
a (Lao) = q+, and kinematic boundary conditions, V = 
V- + dya{0) and V — + dya{L^). 

The ensuing force velocity relations are shown in Fig 10 
for different values of K. Qualitatively, these curves are quite 
similar to their analogs in the mean field model, in particu- 
lar, the negative mobiUty regimes persist for sufficiently small 
K. However, the problem with the divergence of the cell 
length at finite Q does not disappear which means that Kelvin- 
Voigt regularization of tension is weaker than in the mean field 
model. The reason is that the linear stress density dependence 
in the bulk does not penalize sufficiently the infinite stretching 
of the gel layer 

We observe, however, that in the framework of linear 
Kelvin- Voigt elasticity the density singularities shown in Fig|6] 
(b) disappear, which means that this model regularizes infi- 
nite compression adequately, something the mean field model 
could not accomplish. This suggests that the Kelvin- Voigt 
model and the mean field models show complimentary fea- 
tures and should be used in combination. 



p{y)= {dya{y)-V)J,'' 



du 

a„(T(«)-v 



C. Maxwell elasticity 

In contrast to two elastic regularization schemes consid- 
ered above. Maxwell model associates elasticity with fast time 
scales. In the interpretation of this model with co-rotational 
(Jaummann) convective derivative 1281 l52l l53l the dimen- 
sional problem can be written as 

dtp + dx{pv) = 

= d,a (15) 
{ri/E)dtcr + vd^cr + cr = X + ild^v 
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where E is the (infinite frequency) elastic modulus. One 
can see that in this setting the mechanical problem decouples 
again from the mass transport problem. 




Figure 11. Force- velocity curves for the Maxwell model with dif- 
ferent Ai,2,3 ~ {0, 0.1, 1}. Other parameters are po = 1, — 1.7 
and v+ — 2. The minimal model is recovered at A = 0. 



A single dimensionless equation describing steady state 
(TW) regimes takes the form 



\dya{dy(T — V) 



dyya 



(16) 



where the new non dimensional parameter is A = rj/E and 
experimental data suggest that A = 0.02 - 0.2 130115111601 - 
The equation ( [T6| ), which is nonlinear in contrast to what 
we have had in the minimal model, must be again supple- 
mented by two mechanical boundary conditions, (7(0) = (/_ 
and a{Loo) = q+, and two kinematic boundary conditions, 
dya{0) — V = —v^ and dya{Loo) — V = — w+. 



The nonlinear boundary value problem ( 16 1 was studied nu- 
merically and in FigjTTjwe show the force velocity relations 
corresponding to different values of A. One can see that the 
negative motility regimes survive at finite A which suggests 
that the qualitative behavior observed in the minimal model is 
stable under this regularization. 

Introducing Maxwell elasticity, however, fails to regularize 
the infinite stretching singularity. Moreover, Maxwell model 
does not allow for static equilibria describing the rest state of 
a cell and there is a numerical evidence that the problem with 
infinite localization of mass also persists. We can then con- 
clude that at least in the study of steady motility regimes the 
combination of Kelvin- Voigt and mean field elasticity should 
be preferred to the use of Maxwell elastic regularization. 



IV. MISCELLANEOUS 

To check robustness of our predictions we study in this Sec- 
tion three different extensions of the minimal model not deal- 
ing with elasticity. The first extended model allows for inho- 
mogeneous friction, in the second model depolymerization is 
assumed to be taking place everywhere in the bulk of the cell 
body and in the third active contractile pre-stress becomes a 
function of actin density. 



A. Inhomogeneous friction 

Assume that in the steadily moving cell the friction coef- 
ficient f is graded from rear to front. For instance, ^ may be 
viewed as proportional to the steady state density of focal con- 
tacts which are known to concentrate in the frontal part of the 
advancing lamellipodium i30]l63]l64l . 

More specifically, suppose that f = £,k{z) where z = 
y — Lao/ 2 and — Loo/2 < z < Loo/2. The dimensionless 
mecharrical equations describing TW regimes takes the form 



dzQ- 

k{z) 



1. 



(17) 



While this equation is still linear, it now has a variable coef- 
ficient. The mechanical boundary conditions remain the same 
as in the minimal model (t(— Loo/2) = q- and (t(Loo/2) = 
(7+, but the kinematic boundary conditions get modified 



V 



d-.(7 



(-Loo/2) + f_ 



d-.a 



(Loo/2) +?;+. 



A semi-explicit solution of the resulting Sturm-Liouville 
problem can be expressed in terms of two linearly indepen- 
dent functions A{z) and B{z) solving the following elemen- 
tary sub-problems [65J; 



A'(-L, 



A" 

3/2) 



k{z)A 
1 = A'(Lo 



./2), 



and 



S" = k{z)B 
B'(-Loo/2) = 1 and B'(Loo/2) = -1. 



By using the functions A{z) and B{z), we can write the force 
velocity relation in the following explicit form 



AV ^ {2 - Qe)[A] + Q{A} 
V = V^+'-^[B] + ^{B}, 



where 2[f] - /(Loo/2) - /(-Loo/2) and 2{/} - 
/(Loo/2) + /(-Loo/2) . 

Suppose, for instance, that k{z) = 1 + 9ki{z) , where 6 is 
a small parameter and the function ki{z) is odd. Then, in the 
lowest order in 9 we obtain: 

J [A] = = -tanh(Loo/2) 

i {A} ^ -[B] = -^-4^/W2/^sinh(2z)«:i(z)d. 
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The resulting force velocity relation can be written semi- 
explicitly 



AF= (2-Qe)tanh(^)- 



2tanh(- 



L /2 

Observe that if the integral /q °° sinh(2z)Ki(z)(iz in ( 18 i 
is positive which means that if there is a frictional bias at tEe 
front, the cell will have larger length and will move with larger 
velocity than in the minimal model with 6 — 0. If, instead, 
the friction is stronger at the back, the cell will have smaller 
length and will move slower than in the minimal model. These 
results are compatible with the observation that adhesion com- 
plexes predominantly position themselves at the front of the 
moving cell ll30l l63l l64l which can then be interpreted as an 
optimization of velocity. 




Figure 12. Force velocity relations in the case of inhomogeneous 
friction with S as a parameter. The minimal model is recovered at 
6» = 0. 



In Fig 12 we show numerical results for finite values of 9. 
To ensure that the concentration of adhesive complexes at the 
front is four times larger than in the back, which is plausible 
for keratocytes f30l, we must take 9 — 8. From Fig 12 



we 



see that at this level of inhomogeneity the general shape of the 
force velocity curves remains qualitatively the same as in the 
minimal model. 



B. Delocalized depolymerization 



As a result the density distribution does not affect the force 
velocity relation and different models of actin transport can 
be made compatible with the same force velocity relation. 

To exploit this idea we consider in this Section a version of 
the minimal model where depolymerization is not localized at 
the rear front of the cell. Instead, we assume that depolymer- 
ization takes place everywhere in the bulk of the cell as it is 
suggested by observations |66 -69|. 

As a first step we exclude localized depolymerization by 
putting V- = 0. Then we modify the mass conservation equa- 
tion by adding a source term. If we make the simplest as- 
sumption that the rate of depolymerization is a linear function 
of density we obtain 



(19) 



The coefficient 7 can be estimated in the range 0.01 — 0.05s ^ 
Il30l|3l]|70l. 

The proposed reformulation of the minimal problem affects 
the velocity distribution only through the specific choice of 
one of the kinematic fluxes (condition — 0) which gives 
AV = w+ and Vm — v+/'2,. The explicit solution (|3j and 
the general formulas (jSj) remain valid. The stability condition 
takes the form 

Q<v+<2-eQ. 

One can see that the assumption ~ 2, which we used 
throughout the paper to illustrate the results obtained in the 
minimal model, is no longer adequate in the case of pulling. 
We therefore assume a smaller value ~ 1 which is also 
plausible in view of |28 , 31 , 67 , 68 1. In the presence of elas- 
ticity penalizing infinite stretching such re-scaling is not nec- 
essary. 



The treadmilling boundary condition for ( 19 1 can now be 
written in the form 



pil+it),t) 



7M 



(20) 



It ensures that the total mass remains constant 



dM 



-^M + p{l+{t),t)v, 



0. 



We can again absorb M into the scaling of p by using dimen- 
sionless variable p/po- The ensuing nondimensional problem 
depends on the new parameter 7 = 777/x which can be esti- 
mated in the range ~ 0.5 — 2.5. 

The dimensionless equation describing the TW regimes 
takes the form 



dy{p{v-V)) = -jp 



(21) 



where we recall that y — x — Vt. If we now introduce the 
treadmilling mass flux m — —7, we may write the solution of 
( 2 1 1 explicitly 



We recall that in the minimal model the mass transport p^y-j 
problem is fully decoupled from the force balance problem. 



v{y) - V 



exp 



v{u) - V 



du 



(22) 
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Figure 13. Density distribution in the pulling and pushing regimes 
for the model with delocalized depolymerization. Parameters: v+ = 
1 and 7 = 1.5 



Here the pre-exponential factor is exactly the same as in the 
minimal model ( [T3] l while the new exponential term describes 
modulation due to distributed depolymerization. 

To illustrate the role of this term we now show that in this 
new setting we can obtain a peak of density at the front of the 
cell and a decay in the back which is the pattern typically ob- 
served in moving cells lfT8ll67ll68]| . To this end we explicitly 
compute the values of actin density at the front piLaa) and 
at the back p(0). The first of these quantities can be found 
directly 



J_ 



To find the second quantity we need to use an asymptotic de- 



velopment of the integral term in ( 22 1 at small y 



p{y) 



1 



A+ -2A+ sinh(Loo)' 



where A+ and A^ are defined in Q. From this formula we 
see that in the relevant range 7 > A^ (the assumption 7 > 1 
ensures this inequality) we obtain 



p(0) = 0. 

Notice that now we have p(0) < p{Lo 
model we always had p(0) > p{L^). 



while in the minimal 



In Fig 13 we choose 7 = 1.5 and v^^ = 1 and show the 



typical density profiles for both pure pushing and pure pulling 
cases. One can see that even for unloaded cells the distributed 
depolymerization leads to a global decay of actin density from 
the front to the back. Interestingly, this decay may not be 
monotone. It is clear, however, that the fine structure of the 



density at the front of the cell is dominated by microscopic in- 
teractions of the actin network with the membrane and cannot 
be captured by our simplified model. 



C. Density dependent contraction 

The fact that contractile pre-stress depends on both actin 
and myosin densities is well established ll36l |7T1 - I74l . There- 
fore the assumption that the corresponding active term is equal 
to a constant is one of the main weaknesses of the minimal 
model. 

A rather general theory of active gels implying coupling be- 
tween the active prestress and the transport of different com- 
ponents of actomyosin network has been developed in liTSl - 
1221 ■ While it was shown that this coupling induces a rich va- 
riety of dynamic behaviors, the problem was typically studied 
in a fixed domain and without external loading. In particular, 
the issue of the force velocity relation was not addressed. 

In this Section we study force velocity relations under the 
assumption that contractile stresses depends on cytoskeletal 
density representing actin filaments with the same orientation. 
At the same time we neglect the important coupling of active 
stress with motor density studied in 1321 l36l l74l . 

We begin by writing the system of coupled equations of the 
model where it is convenient to distinguishing three subprob- 
lems: 

Force balance: 



-dxxO- + cr = aaip) 

(j{l-{t), t) = (7_ and a{l+{t),t) = q+ 
Polymerization/ depolymerization : 



(23) 



i+ = «H 

L(0) = 



+ dx<j{i+{t),t) 
f < l+{0) = /o 



(24) 



Actin transport: 



dtP + dxipdxcr) — 
p{l^{t),t)v- = p{l+{t),t)v+ 
pix,0) = pi{x). 



(25) 



The force velocity relation can be obtained by solving these 
equations numerically for different loading conditions and 
tracing the solutions till they approach different traveling 
wave regimes. 

In our numerical experiments we observed that by assum- 
ing a linear dependence of active prestress on actin density 
o'a(p) we do not reach TW regimes and obtain instead os- 
cillatory modes of cell motility 177 j. The situation changes 
if we assume that this dependence is non-monotone which 
agrees with microscopically motivated models considered in 
II7TI473I l78l ■ According to these models at small actin densi- 
ties more filaments allow more motors to bound and to induce 
contractile stresses, however, there is a density threshold after 
which compaction of the network prevents further increase of 
the contractile stresses. By accepting this reasoning we used 
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- C=1 ,p^=4 

- C=1 , p^=ra 

- C=4, p =2 



b 



Figure 14. Contractile stress Oa as a function of actin density p for 
different choices of parameters ps and C taken from L78J . Minimal 
model is recovered at C = 1, ps = oo. 



the function (Ta{p) proposed in ll78l which in dimensionless 
form can be written as 



(26) 



Here C is a constant and ps is the actin saturation density, see 

Figjg 

In our numerical experiments the initial location of the cell 
boundaries was at = 0, = 1. The initial density distri- 
bution was chosen to be 



1 + u_ / w+ 



1 



v+ 



1 X 



which is the simplest way to satisfy the boundary conditions 

1° 

and the requirement that p — 1. By varying the initial 

data we could reach different traveling wave regimes and in 
this way recover the full force-velocity relation. 



The results are presented in Fig 15 As we see, the imposed 
coupling does not destroy the fundamental difference in con- 
vexity properties between the force velocity curves in pushing 
and pulling regimes. We also observe that in the parametric 



regimes presented in Fig 15 the steady density profiles always 
lie on the decreasing limb of the density-contractility curve 
shown in Fig 14 Further studies are needed to understand 



this phenomenon as well as other effects including oscillatory 
and stick slip type non equilibrium steady states li75-,-,-77J . 



V. THE EFFICIENCY OF CARGO TRANSPORT 

The simplicity of the minimal model allows one not only 
to obtain explicit force velocity relations but also to study the 
energetics of a self propelling cell carrying a cargo. 

If we multiply the force balance equation ([T]i by u(a;, and 
use the constitutive realtion (|2| we obtain the global energy 
balance equation 



;+(t) 



;+(t) 



;+(t) 



i'+(t) 



In this equation we can identify the following terms 

1. Df = //^(j*"* v"^ > 0, dissipation rate associated with 
surface friction. 




;:ES 




0.05 1 15 2 25 0.3 35 



Q 



Figure 15. Typical force velocity relations for the model with actin 
dependant contractility. Parameters are u_ — 1.7, w+ = 2. Inserts 
show the steady states density profiles located on the decreasing limb 
of the density-contractility curve, see Fig[T4] Minimal model is re- 
covered at C = 1, Ps = oo. 



2. Dy = JI^^^^ (dxv)'^ > 0, dissipation rate associated 
with bulk viscosity. 

3. Pc = — Ji^(t) > 0, rate of energy consumption by 
the contractile mechanism. 

4. Pp ~ {qv)- — iqv)+ > 0, rate of energy consumption 
by the protrusion mechanism. 

5. A — (ql)- — {ql)+, the power expanded against the 
external forces. 

In the case of TW regimes all these terms can be computed 
explicitly. In particular, by using nondimensional quantities 
we obtain, 

P, = AV, Pp = QV,n - A = QV. 

In these notations we can write the energy balance in the form 

Pp + Pa = A + D, (27) 

where D — Df + Dy. The mechanical efficiency of cargo 
transportation can then be defined as follows: 



A = W/H. 



(28) 



Here the numerator W describes the useful work and may in 
addition to A contain an additive Stokes term Pg = LV^, 
which is nonzero even in the absence of the cargo ||79]48T1 . 
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However in our problem this correction can be shown to be 
small and will be neglected. 

The denominator H > describes external energy supply 
associated with ATP hydrolysis which drives the motility pro- 
cess. It is clear that H must include the power Ha — Pp + Pc 
exerted by active forces on the constraining environment. It 
should also contain the 'maintenance' term Hjj which ac- 
counts for energy consumption required to sustain the active 
state in the absence of macroscopic motion ll82l . By using ter- 
minology introduced in |28, 83 1 for weakly non-equilibrium 
regimes we can identify Ha and Ha with the terms that are, 
respectively, linear and quadratic in the measure of chemi- 
cal non-equilibrium A/i. In what follows, we neglect the 
'quadratic' term H]j (dealing with degrees of freedom that are 
invisible in our macro-scopic model) comparing to the 'linear' 
term Ha- 
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Figure 16. Rate of energy consumption and efficiency as functions 
of the load in the cases of pure pulling and pushing loading modes. 
The corresponding force velocity relation is shown in Fig|2] Driving 
parameters are v- — 1.7 and v+ = 2. 



displays a usual single maximum, in the case of pulling the 
efficiency-load relation becomes bi-modal. The two maxima 
can be identified with protrusion dominated and contraction 
dominated motility mechanisms. Such bi-modality may carry 
biological advantages allowing a the cell to switch back and 
forth between two highly efficient regimes by controlling, for 
instance, the friction experienced by the nucleus. 



A 0.4 




Figure 17. Efficiency as a function of the load in the elasticity- 
regularized model in pure pushing and pulling regimes with different 
fci,2,3,4 = {0,0.01,0.1, 1} and Lo = 1. Experimental data suggest 
that fe = 1 - 10 (see (23l|24l). Parameters: = 1.7 and v+ = 
2. The corresponding force velocity relations are shown in Fig[8] 
Minimal model is recovered at A; = 0. 



If we augment the minimal model by adding mean field 
elasticity (see Section [lll| the energy balance equation takes 
the form 



P„ + Pc+Pe=A + D, 



By using these simplifying assumptions we can write 



A 



A 



Pn + Pr 



(29) 



By using the fact that D > one can showthat < | A| < 1. 



In the TW limit the efficiency ( 29 1 can be computed explicitly 



A = 



QV 



AV + QVrn - eQAV/2 ' 



In Fig 16 we show the efficiency A and the energy con- 
sumption rate Pp + Pc as functions of the total load Q for 
the TW regimes presented in Fig|5] First of all we observe 
that the divergence of the cell length in the minimal model 
in the pulling regimes with Q > Qc does not lead to singu- 
lar behavior of any of the energetic measures. Second, we 
notice, that while in the case of pushing the function A{Q) 



where the new term P^ describes the power exerted by protru- 
sion mechanism against the elastic 'spring' 

^0 

This term remains nonzero in TW regimes where L = and 
it should be added to the denominator in the expression of the 
efficiency 



A 



A 



Pv + Pc+ Pe 



(30) 



In Fig 17 we show the k dependence of the efficiency ( [30| . 
One can see that the two peak structure of the function A{Q) 
survives in the regularized model till a threshold in k, signi- 
fying also the disappearance of the negative mobility range, is 
reached. 
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VI. ALTERNATIVE DRIVING MODES 

In the minimal model we used an assumption that the pro- 
cess of cell motility is driven by the kinematic fluxes charac- 
terized by parameters v+ and v- . This assumption, illustrated 
in Fig. 
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means that we impose separately the velocities of 
polymerizing (arriving) and depolymerizing (departing) mass 
points, see also I71|20l[l8l[30l^. The fact that nothing is said 
about the densities of the arriving or departing material allows 
one to decouple the mechanical problem from the mass trans- 
portation problem and makes the analysis fully transparent. 
This transparency, however, comes at a cost. 



-wmwvMu- 



i-{t) 



Figure 18. Schematic structure of the treadmilling cycle showing 
different densities of arriving (polymerizing) and departing (depoly- 
merizing) material. 



First, it is clear that the treadmilling is characterized by only 
one parameter, the mass flux m, so by fixing two parameters 
Vm and AV we are implicitly constraining both treadmilling 
and contraction. This is also clear from the fact that parameter 
AV = w+ — u_ serves as a measure of energy consumption in 
the contraction mechanism. Second, by prescribing the kine- 
matic fluxes v+ and u_ we have no direct control of the tread- 
milling mass flux. As a result we encounter singular regimes 
with rn — Q which leads to either infinite mass localization 
inside the cell or to infinite spreading of the cell body. Third, 
by focussing on kinematic fluxes we do not put any restric- 
tions on the energy consumption required to sustain different 
active mechanisms which appears to be a natural biological 
constraint. 

Notice also that the problem setting where driving is per- 
formed through parameters v+ and V- contains an implicit 
assumption that the material arrives with a particular density 
(particular structural organization). Another implicit assump- 
tion is that the departing material has a density which depends 
on the activity of the contractile machinery. While these as- 
sumptions are plausible, they may not be the most natural ones 
from the biological point of view. Even more importantly, as 
we have shown in the previous Sections, these assumptions 
necessarily lead to singularities. 

In this Section, to further challenge the robustness of 
our conclusions about the asymmetry between pushing and 
pulling, we consider an alternative modality of driving by im- 
posing constraints on energetic rather than kinematic parame- 
ters. The main difficulty in dealing with non kinematic driv- 
ing schemes is that they couple the mechanical and the mass 
transport problems already in the minimal setting. 
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Figure 19. (a) Various pulling regimes in the parameter space 
(m, Pc). The singular regimes correspond to the lines ZX (infinite 
localization) and to the point X(R) (infinite spreading), (b) Density 
localization along the path indicated in (a) by the solid line which 
ends with the formation of a singularity at point W. The loading is 
e = 1 and Q = 1.6. Analogue of Fig|6] 



More specifically, we assume that the cell controls the 
treadmilling rate, characterized by the total mass flux rh < 0, 
and the energetics of the contraction process, characterized by 
the consumed power Pc = AV. The advantage of this new 
parametrization is that protrusion and contraction can now be 
controlled independently. If we choose the pair (m, Pc) as the 
parameters instead of (Vm, AV), we again obtain stable TV 
solutions given that Pc < 2 — Qe and m < 0. The proposed 
driving mode is in fact equivalent to the kinematic driving 
mode in the TW regimes because the Jacobian of the trans- 
formation — ^ {Pc{{v-,v+)),m{v^,v+)) 




Jo 



{v{y)-Vy 



ilo 



dy 
v(y)-V ' 



> 



1 



> 



is strictly positive for < Loo < oo. 

By using the parametrization (to, Pc) we can easily avoid 
the density localization phenomenon illustrated in Fig|6]with- 
out introducing elasticity. To illustrate this point we show in 
Fig 19 the pulling TW regimes in the parameter plane (to, Pc) 
where we again distinguish between regimes where cell car- 
ries the cargo (V > 0) and regimes where it is dragged by the 
cargo (V < 0). Fig 19 has to be compared with Fig|6]where 
the same regimes are shown in the (Vm, AV) space; the only 
difference is that now the line XR corresponding to regimes 
with Loo ~ oo collapses on a single point X{R). It is clear 
that if the treadmilling flux is prescribed so that to 7^ 0, the 
singularities associated with the line ZX{R) in Figj6]are au- 
tomatically excluded. 



In Fig. 20 we show the force velocity relations in the mini- 
mal model with prescribed {rn, Pc). One can see that the qual- 
itative difference between pushing and pulling endures in this 
new setting, moreover, we again observe regimes with neg- 
ative mobility. It is interesting that by fixing parameters Pc 
and 771 we induce a dependence of the polymerization and de- 
polymerization rates (w_ , w+) on Q (see the inserts in Fig. |20ll 



which agrees qualitatively with the trends suggested in ED i 



based on the polymerization ratchet model. We also note that 
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Figure 20. Force velocity relations in pure pushing and pulling 
TW regimes when driving is performed by imposing rh — —6.1 and 
Pc = 0.3. Inserts show the ensuing dependances of v+ and V- on 



at sufficiently strong pulling loads Q > — {2 ~ P^) /e, the 
cell length L^o diverges which suggests that also in the case 
of non-kinematic driving the minimal model should still be 
elastically regularized. 

Finally we remark that instead of the pair (m, Pc) we 
could also prescribe another set of energy related parameters, 
for instance, {Pp,Pc). Ultimately, the choice of the driving 
mode requires an understanding of the microscopic side of 
the model and the answer may depend on the type of the cell, 
the environment and the regime of loading. 



VII. CONCLUSIONS 

In this paper we used the simplest model of a crawling cell 
to study an interplay between contraction and protrusion re- 
quired to sustain and carry various cargoes. The model de- 
scribes a layer of active gel subjected to external forces. It 
extends previous studies focussed predominantly on the be- 
havior of unloaded active media or on problems with fixed 
boundaries. 

By using an analytically transparent framework provided 
by the minimal model we demonstrated for the first time that 
contraction and protrusion mechanisms can interchange their 
roles as one varies the dipole component of the external load. 
Our model predicts a possibility of a relatively sharp transi- 
tion between protrusion dominated motility and contraction 
dominated motility in response to an increase of the pulling 
force. This transition has a macroscopic signature and can be 
in principle identified experimentally with a negative mobility 
range on a force-velocity curve. 

The advantage of the minimal setting is that it delivers ex- 



plicit steady state solutions describing asymmetrically loaded 
self-propelling cells and allows analytical access to their sta- 
bility. Only in such prototypical framework the competition 
between contraction and treadmilling can be studied in the 
transparent form without any geometric effects obscuring the 
interplay between competing active ingredients of the model. 
The simplicity of the model allowed us to elucidate active 
adjustment of the force producing machinery to the subtle 
changes in the character of external loading. The possibil- 
ity of such adjustment implies that 'pushers' both collaborate 
and compete with 'pullers'. 

The augmentation of the active gel model involving elastic 
stresses in addition to viscous and active stresses was found 
to be essential for the removal of singularities inherent in the 
minimal model. In particular, mean field elasticity appears to 
be the most universal way of introducing a resting configura- 
tion even when treadmilling is absent and to deal with infinite 
spreading at finite pulling loads. Instead, Kelvin- Voigt visco- 
elastic model while also ensuring existence of static configu- 
rations and removing infinite density localization, fails to se- 
cure the finite length of the cell in the whole interval of applied 
pulling loads. We conjecture that a combination of mean field 
and Kelvin- Voigt elastic terms in the system of equations de- 
scribing active gels is sufficient to fully regularize the minimal 
model. 

To make definitive predictions about the feasibility of the 
negative mobility regimes, focussed measurements of the ef- 
fective stiffness associated with different elastic structures of 
the cell are necessary. In the situation when elastic coupling 
strongly affects the force velocity relations, studying kinetic 
relations for differently loaded cells may be the way to furnish 
a set of independent bounds on such stiffness. It should be em- 
phasized, however, that our conclusions regarding convexity- 
concavity structure of the force velocity relation are much less 
sensitive to the value of the stiffness than the very existence of 
the negative motility regime. Thus, pushing and pulling force 
velocity relations remain qualitatively different even when the 
negative motility regime disappears. 

Perhaps our most intriguing finding is that the fine structure 
of the force-velocity relation may depend on the modality of 
external driving and we argued that kinematic driving may not 
be the only physically and biologically natural choice. In par- 
ticular, we suggested that instead of the rates of polymeriza- 
tion and depolymerization, the cell may be controlling the en- 
ergy supplies required for the functioning of contraction and 
protrusion mechanisms. We have shown, however, that while 
the detailed shape of the force velocity relation depends on the 
choice of the driving mode, its loading-sensitive convexity- 
concavity structure is a robust feature of the model. 
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